APS/123-QED 



Lattice Boltzmann - Langevin simulations of binary mixtures 
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We report a hybrid numerical method for the solution of the model H fluctuating hydrodynamic 
equations for binary mixtures. The momentum conservation equations with Landau-Lifshitz stresses 
are solved using the fluctuating lattice Boltzmann equation while the order parameter conservation 
equation with Langevin fluxes are solved using the stochastic method of lines. Two methods, 
based on finite difference and finite volume, are proposed for spatial discretisation of the order 
parameter equation. Special care is taken to ensure that the fluctuation-dissipation theorem is 
maintained at the lattice level in both cases. The methods are benchmarked by comparing static 
and dynamic correlations and excellent agreement is found between analytical and numerical results. 
The Galilean invariance of the model is tested and found to be satisfactory. Thermally induced 
capillary fluctuations of the interface are captured accurately, indicating that the model can be used 
to study nonlinear fluctuations. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

Thermal fluctuations are an essential part of the 
physics at mesoscopic length scales in fluid mechani- 
cal problems. For instance, thermal fluctuations pro- 
duce Brownian motion in colloidal suspensions, confor- 
mational fluctuations of polymers and membranes, cap- 
illary waves at fluctuating interfaces, and critical opales- 
cence in binary mixtures. A consistent mesoscopic de- 
scription of such phenomena follows from the equations 
of fluctuating hydrodynamics. The first instance of such 
a description was the fluctuating Navier- Stokes equa- 
tions of Landau and Lifshitz [l| . Similar equations were 
then introduced to study the dynamics of order param- 
eter fluctuations in criticalphenomena, as reviewed by 
Halperin and Hohenberg 0]- The coupled fluctuating 
equations of motion for the momentum and order pa- 
rameter are known as model H in their classification. 

The model H equations describe the fiuctuating hy- 
drodynamics of a conserved order parameter and the 
conserved momentum density g = pu, where p and u are 
the total density and the local fiuid velocity. To ensure 
conservation of local densities, fiuctuations are incorpo- 
rated as random stresses in the momentum equation [l| 
and as random fluxes in the order parameter equation [3| . 
At equilibrium, these random fluxes are constrained by 
fluctuation-dissipation theorems, which relate their vari- 
ances to the kinetic coefficients in the equations of mo- 
tion. The fiuctuation-dissipation theorem (FDT) ensures 
that the dynamical equations give rise to a Gibbs distri- 
bution for the fiuctuating variables, as required by equi- 
librium statistical mechanics. Thus, together with the 
conservation laws, the FDT is an important constraint 
in the model H equations. 

The model H consists of non-linear stochastic partial 
differential equations which admit no analytical solu- 
tions, requiring, therefore, numerical methods of solu- 



tion. Numerical methods which proceed by discretising 
the equations of motion on a lattice must ensure, at least, 
that the conservation laws and the FDTs are obeyed. 
This requires care as naive discretisations often violate 
the FDT, leaving degrees of freedom incompletely equili- 
brated, and therefore, without a Gibbs distribution [j-Q . 

In this paper, we solve the model H equations by 
combining the fluctuating lattice Boltzmann equation 

SLBE) 0, H] with a stochastic method of lines (SMOL) 
[loj using both finite difference and flnite volume dis- 
cretisations [lT| . The formulation ensures conservation 
of local densities to machine precision, and a correct bal- 
ance between fluctuation and dissipation for all the de- 
grees of freedom on the lattice. We expect our method 
to be widely applicable to problems in binary mixtures 
and other physical systems where model H is applica- 
ble, when thermal fluctuations form an essential part of 
the physics [l^ - [T5| . Hybrid methods have been devel- 
oped in the literature in different contexts, for example 
in case of dynamics of binary complex fluids [l^ , but 
without considering thermal fluctuations. We deal with 
fluctuating hydrodynamics of binary fluids in detail here. 
Alternative schemes based on finite volume methods have 
also been used to simulate the fluctuating hydrodynam- 
ics of single component fluids [ist and reaction-diffusion 
systems [l9|. However, the methodology outlined here 
carries the advantages of the lattice Boltzmann method 
[20| and can be generalized to other problems in fluc- 
tuating complex fluids, for instance, to the dynamics of 
microemulsions [2l| and liquid crystals. 

The rest of the paper is organized as follows. In the fol- 
lowing section we provide a detailed description of model 
H. We review the current understanding of solving of 
these equations in section IIIII In section IIVI - IVl we 
present the numerical method followed by the validation 
results in section IVll We compare our method with pre- 
vious approaches and end with a summary of our work 
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in section rVlII 



II. FLUCTUATING HYDRODYNAMICS OF A 
BINARY FLUID MIXTURE 

We consider a coarse grained model for an isothermal 
binary fluid system, consisting of species / and // with 
local densities ni and nn. The mixture as a whole has 
density p = nj + n//. The order parameter -0, which 
quantifies the local composition, is taken as the normal- 
ized density difference. 



ni - nil 
ni + nil 



(1) 



A. Landau-Ginzburg theory 



the order parameter /i = SF/5ip = Aijj + Bijj^ — KV'^tp. 
The three parameters A, B, and K control the interfa- 
cial thickness and interfacial energy of the mixture and 
after suitable non-dimensionalisations, allow for compar- 
isons with real fluids. The additional stress due to the 
presence of order parameter gradients follows from the 
relation i/^V/i — V ■ cr"^ [1^, which is solved by 



df 



SafS 



K 



(7) 



This additional stress includes the Laplace and 
Marangoni stresses due to a fluid-fluid interface. The 
form of this stress tensor can be motivated on the ba- 
sis of an electrostatic analogy or derived directly from 
Poisson brackets 123 • 



The equilibrium thermodynamics of the fluid is de- 
scribed by the Landau free-energy functional [l^, 



(2) 



Here, ■0 is allowed to vary beyond the limits of ±1 that 
follow from its definition. This "softening" of the or- 
der parameter has no consequence in the thermodynamic 
limit jl^l- The first term represents the local free energy 
density of the bulk fluid, and is approximated as 



(3) 



with A < Q and B > Q. The second term of Eq. [2] in- 
volving the square gradient gives a free energy cost to any 
variation in the order parameter, and is related to the in- 
terfacial tension between the two fluid phases [1^ . Min- 
imization of Eq. 13] with respect to the order parameter 
gives two uniform solutions ip = ± ^ A/ B, corresponding 
to two equilibrium fluid phases. These two phases can 
coexist through a fluid interface. For a planar interface, 
the profile joining the two bulk phases reads 

/'a z 
- tanh - (4) 

where z is the co-ordinate normal to the interface while 



(5) 



determines the interfacial thickness. The excess energy 
associated to this profile with respect to the bulk energy 
provides the interfacial tension 



2 2KA^ 
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B^ 



(6) 



The corresponding chemical potential is given by the 
variational derivative of the free energy with respect to 



B. Model H 

Model H of Halperin and Hohenberg describes the 
coupled dynamics of a conserved scalar order parame- 
ter ip and the conserved momentum density g. The or- 
der parameter is described by a fluctuating Cahn-Hilliard 
equation, known as model B, which includes advection by 
fluid flow, relaxation due to chemical potential gradients, 
and spontaneous thermal fiuctuations. 



(8) 



The mobility M is the constant of proportionality in the 
linear phenomenological law relating the thermodynamic 
fiux of tp to the thermodynamic force V/i. We consider 
M to be a constant, though such an assumption is not 
necessary. Thermal fiuctuations associated with ip are 
introduced through the random flux ^. 

The order parameter dynamics is coupled to a fluc- 
tuating Navier-Stokes equation [l| with additional stress 
densities arising from the order parameter. For a com- 
pressible fluid, the dynamics is governed by 



9tg + V- (ug) = -Vp + vyV^u 



d-2 



d 



-Tl + Vb 

V-0- 



V(V ■ u) 

(9) 



together with the continuity equation for the density. In 
the above, p stands for the isotropic contribution of the 
pressure. & is the random stress introduced by Landau 
and Lifshitz, ■i/'V/i is the order parameter stress rj and 
r]h are the shear and bulk viscosities respectively and d 
is the dimensionality of the system. Qualitatively, these 
equations describe the coupled dynamics of order param- 
eter and flow : inhomogeneities in the order parameter 
generate chemical potential gradients, which in turn pro- 
duce stresses in the fluid. These stresses are relaxed by 
fiuid fiow, which in turn advects the order parameter to 
produce inhomogeneities. 
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The random flux is a zero-mean Gaussian white noise 
whose variance is fixed by the FDT to be 

i,ir,t)ifi{r',t')) =2kTMd^pS{r~r')S{t^t') (10) 

for an isothermal fiuid at temperature T, where k stands 
for the Boltzmann constant. Similarly, the random stress 
is a zero-mean Gaussian white noise whose variance is 
fixed by the FDT to be 

{ac.p{r,t)a^s{r' ,t')) = 2kTr]^p^sS{r - r')S{t - t'). (11) 

where rj^p-ys is the tensor of viscosities formed out of the 
isotropic tensor 5ap and the shear viscosity, 77, and bulk 
viscosity, 77;,: 



= V (Sa-fSpS + SasSpj) + yrjb - -rjj SafjSyS (12) 

For simplicity we assume the same viscosity for the two 
fluid phases. 

In the next section we briefly review previous algo- 
rithms to numerically solve these coupled equations and 
point out why they lead to an incomplete equilibriation 
of both the order parameter and momentum degrees of 
freedom. This drawback imposes severe restrictions in 
the applicability of these algorithms to situations where 
a complete equilibriation is required, a gap which our 
work attempts to flU. 



III. DISCRETISATION AND FDT VIOLATION 

There are ample instances in the literature where a 
naive discretisation of both the momentum Q and order 
parameter 0, 0, HI] equations have led to FDT violations 
on the lattice. An important question, then, is how best 
FDTs. derived in the continuum with respect to appro- 
priate conservation laws, can be implemented in discrete 
space and time. In this section, we present a very brief 
survey of previous numerical schemes, to clarify when 
naive discretisations lead to FDT violations. 

In order to gain insight into the inconsistencies associ- 
ated with the order parameter discretisation, let us con- 
sider a low-order discrete representations of the diver- 
gence of a vector ^ and the Laplacian of a scalar tp, 



(r) = ^t^iCt -iir + c^) 



[V^tA] (r) =^S,^(r-f c). 



(13) 
(14) 



Here, tUi and Qi are weight factors which depend on the 
stencil, i refers to the number of neighboring grid points 
considered, {c^} corresponds to a lattice vector and hence 
r -I- Ci represents the points of the chosen stencil. In a 
Fourier representation, they become 



(q) - ^^.ce^'i-' • |(q) = r(q) • |(q) (15) 



where r(q) and i(q) are the Fourier representations of 
the divergence and Laplacian operators, respectively. It 
is easy to see that r(q) — !• iq and L{q) — > —q^ as 
q —7- for any admissible choice of stencil. In that limit, 
we recover the lattice analogue of the familiar relation 
between the gradient and Laplacian operators, so that 
i(q) = r(q) • r(q). At high wavenumbers, however, this 
relation is no longer true. Indeed, it is violated by all 
standard nearest neighbour stencils [1^ . 

To see how this affects discretisations of the fluctuat- 
ing Cahn-Hilliard equation, we linearize Eq. |S] about a 
state of zero flow, for completely local and harmonic free 
energy (B = 0, if = in Eq. [U, with a mobflity that 
is independent of the order parameter. Discretising and 
Fourier transforming, we obtain 

dMq) = A/i(q)A^(q) + r(q) • ^(q). (17) 

It is evident from Eq. [T7] that fluctuations in the order 
parameter equation will satisfy the FDT of Eq. [TU]on the 
lattice if and only if L{q) = r(q) • r(q). Equivalently, 
the discrete operators should satisfy = V.V in real 
space. Since this is not true for the standard choices 
of the previous operators [2§|, resulting discretisations 
violate FDT. 

To verify the above analysis, we perform simulations 
using the method proposed by Petschek and Metiu Q 
and used, for example, in Q and [28| . Their method is 
essentially the one outlined above, with specific choices 
of the gradient and Laplacian. Simulations are carried 
out on a 32 X 32 X 32 domain with a cubic grid and unit 
spacing. Ax = 1 and unit time step At = 1 using a 
stochastic Runge-Kutta algorithm [301 ■ We compare the 
theoretical value of the Fourier mode amplitudes of order 
parameter as given by the Gibbs distribution 



Mq)l^) = ^ 



(18) 



[V^T^] (q) = J2 ^^e^''"^' V(q) = i(q)V'(q). (16) 



with our simulation data. We define the equilibrium ra- 
tio ER as the ratio of simulated values to the theoretical 
value. If all Fourier modes are in equilibrium, the ER will 
be unity as dictated by Eq. [THl The results obtained are 
displayed in Fig. [T] As can be seen, the difference from 
the expected theoretical value of ER = 1 is quite signif- 
icant : the match is restricted to only small wave num- 
bers and clearly shows the breakdown of FDT at high 
wavenumbers. Having identified the spatial discretisa- 
tion as the main source of error in FDT violation on the 
lattice, we will analyze in the next section how to circum- 
vent it for a scalar order parameter. Its generalization to 
vector and tensor order parameters is straightforward. 

Fluctuations have also been included in lattice Boltz- 
mann equation (LBE) to recover fluctuating Navier- 
Stokes equations. Ladd Q proposed a modiflcation of 
the LBE with the addition of fluctuating stresses. A 
Langevin interpretation of the Boltzmann equation then 
yields the equations of fluctuating hydrodynamics [l|. 
However, as was pointed out in [7|, Ladd's method en- 
sures thermalisation only in the small wave number limit. 
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space is discretized into n components in d dimensional 
space, the discrete form of the fluctuating Boltzmann 
equation reads HI] 

dth + c. • V/, + [F • Ve/]. = -Y, L,, if, - fj) + C. (19) 



where F(x, t) is an effective force density, Ci(x, stands 
for the fluctuations in the populations, and is the 
discrete form of the collision integral and is related to 
the fluid viscosity. The moments of the single particle 
distribution function fi, defined at lattice node x with 
velocity at time give the fluid mass, momentum and 
stress densities: 



FIG. 1. (Color Online) Equilibrium ratio (ER) according to 
Eq. [T8]as a function of wave vector magnitude q alone the di- 
agonal X = y = z from using a conventional method [a] based 
on finite difi'erence discretisation for both the divergence and 
Laplacian operators. Simulation results show significant dif- 
ferences with theoretical predictions at large wavenumbers. 



This was resolved by relating thermal fluctuations to all 
sources of dissipation associated with the collision oper- 
ator in the lattice Boltzmann equation, leading to ther- 
mal equilibrium for all modes, including the ghost modes 
31 1 . This was confirmed subsequently in M. The fluctu- 
ating lattice Boltzmann equation (FLBE) p, [1] provides 
a consistent lattice discretisation for the Navier-Stokes 
equations and is the approach we shall use in this work. 

The FLBE approach has recently been generalized to 
hydrodynamic fiuctuations of non-ideal gases [13, H^, but 
only a few studies have addressed thermal fluctuations 
in binary mixtures in the context of LBE. Noise driven 
spinodal decomposition was studied in by combining 
Ladd's fluctuating LBE with a fluctuating kinetic equa- 
tion for the order parameter. However, this method does 
not respect FDT for either the momentum or the order 
parameter. Extending this to binary fluids maintaining 
FDT is considerably more difficult and so we prefer the 
alternative hybrid method described below. 



(20) 



^Cifj — c^Sa/s- The collision operator 



ij controls the relaxation of fj to equilibrium, /j*. We 



where Q^ap 
L 

can take advantage of the hyperbolic character of the 
FLBE and use the method of characteristics to evolve 
Eq. [19] over a finite time step. When accounting for the 
effect of forces and fluctuations in the evolution of fi, 
it is convenient to introduce the auxiliary distribution 
function 



/,(x,t) = /i(x,t) 



At 



i?»(x,t) 



(21) 



in terms of i?i(x, t) = ~ J^j ^ijifj — fj) + which rep- 
resents the effects of collision, forcing and thermal fluc- 
tuations, see Appendix |B] For a single-time relaxation 
operator, Ly = Sij/r, the hydrodynamic variables are 
related to the auxiliary distributions as 



4=0 

n 

pVa = E ~^ 

i=0 
n 

Sap = E/ fiQiap 



At 

At/2 



(22) 
(23) 



4=0 



1=0 



IV. FLUCTUATING NAVIER-STOKES SOLVER 



■ + At/2 

+ pVaVfi + T{VaFj3 + FaVfj) + T E OQia;? 



(24) 



We use the FLBE method for solving the fluctuating 
Navier-Stokcs equations. The FLBE introduced in [3| 
needs to be modifled to include force densities, which in 
the hybrid method, are the divergences of order param- 
eter stresses. Since the combination of noise and exter- 
nal force densities [IJl modifies the moment relations be- 
tween the distribution functions and the hydrodynamic 
variables, we discuss now the main new features of FLBE 
and provide a detailed and self-contained derivation in 
Appendix [B] 

In a standard DdQn LBE model where the velocity 



where the equilibrium distribution, can be recon- 
structed from p and pw. In Eq. ([24]) X)i=o dQiap is 
the fluctuating contribution to the stress. 

The effective force density is the divergence of the order 
parameter stress 



F = V • cr"^ = VV^ 



(25) 



which can be verified using Eq. [7] To compute this 
force density we use a symmetrized, second order ac- 
curate nearest-neighbor central difference stencil for the 
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gradient 

V/i(a;, y, z) = i [fi{x + 1, y, z) - fi{x ~ 1, y, z)] x 

+ ^ [^(a;, y + l,z)- fi{x, y -I, z)] y 

+ \ [y^ix, y,z + l) - n{x, y,z- 1)] z (26) 

and the Shinozaki-Oono discretisation [s^ of the Lapla- 
cian (Eq. IA4p to calculate V^t/i in the chemical potential. 

V. FLUCTUATING CAHN-HILLIARD SOLVER 

We use a stochastic method of lines (SMOL) discreti- 
sation to solve the fluctuating Cahn-Hilliard equa- 
tion for the order parameter. Since it does not contain a 
pressure term which acts as a Lagrange multiplier in the 
incompressible Navier-Stokes equations, there is no par- 
ticular benefit in using a kinetic algorithm with its large 
number of degrees of freedom in solving for a single scalar 
variable. Here, we adopt a semi-discretisation strategy 
J discretising the spatial variables to obtain a set 
of coupled stochastic ordinary differential equations. The 
spatial discretisations we propose ensure that the conser- 
vation law is respected to machine precision and that the 
fluctuation and dissipation are in balance for all wave vec- 
tors. We propose a finite-difference and finite- volume dis- 
cretisations, discussing their relative merits below. The 
temporal integration of the resulting stochastic differen- 
tial equations is done using a Runge-Kutta algorithm 
proposed recently by Wilkie [s^l- This is a straightfor- 
ward generalization of the deterministic Runge-Kutta al- 
gorithm where the noise is held constant through the in- 
tegration step. The methodology may be improved using 
implicit schemes to increase the accuracy. 

A. Finite difference method 

To proceed towards a discretisation of the fluctuating 
Cahn-Hilliard equation which preserves FDT, we write 
the order parameter evolution equation in Fourier space 

dt^iq) = ML{q)Jl{q) + T{q) ■ ^(q) (27) 

assuming a constant mobility. Defining the divergence of 
the noise in Fourier space as 77 (q) = r(q) • ^(q), we see 
that it must satisfy 

(77(q, mq',t')) = -2kTMLiq)Siq + q')<5(t - t'). (28) 

Instead of constructing a divergence operator r(q) which 
satisfies r(q) • r(q) = L{q) we directly use the above re- 
lationship to construct the noise in Fourier space. This is 
then inverse-transformed to real space to provide a noise 
which has correlations compatible with the discretisation 
of the Laplacian and the same Laplacian stencil is used 




FIG. 2. Illustration of the stencil used for the numerical tests 
in the finite volume method for a two dimensional case. This 
stencil corresponds to the D2Q9 lattice Boltzmann model. 
Physical quantities ,e.g. -0, V/i, and u, are defined at node 
r which has its neighbors at r-f c^. All fluxes j',^* (diffusive, 
convective and random) are deflned at the mid point of the 
links (r 4- ici) connecting r and r -I- Ci . (See Eq. 1301 - I33|l 

to calculate V^/i and V^V- The generation of noise in 
Fourier space has been used earlier in spectral methods 
[36j to respect FDT in discrete space. 

It is important to ensure as isotropic a discretisation of 
the Laplacian as possible, to avoid artifacts like spurious 
pinning of interfaces by the lattice. We have compared 
in Appendix \X\ four standard finite-difference stencils re- 
ported in the literature, see Fig. [T2]in Appendix \X\ where 
expressions for their Fourier transforms L(q) are also pro- 
vided. The Laplacian of Shinozaki and Oono [S^ is the 
most isotropic one and we use it for our discretisation. 
The advective flux, V- {uip), is discretized using a second 
order accurate, conservative, central difference scheme 

[V-iuiP)] {x,y,z) 

= 2 {["^■^■'K^ + 1'^'^) ~ [uxi^K^ - ^,y,z)} 

+ ^ {[uyi']{x,y + l,z) - [uy^]{x,y - 1, z)} 
+ ^ {[uzip]{x, y,z + l) - [u;,il;]{x, y,z- 1)|29) 

B. Finite volume method 

It is possible to formulate an alternative discretisa- 
tion for the fluctuating Cahn-Hilliard equation, based 
on a finite-volume formulation. Such an approach, us- 
ing fluxes defined on lattice links, has been proposed to 
study the electrokinetic equations in the absence of fluc- 
tuations in [13, • Alternative finite volume schemes 
may also be found in the context of reaction-diffusion 
systems [l^. Specifically, we choose a DdQn cubic lat- 
tice and a set of link vectors {ci} as done usually with 
lattice Boltzmann models. Thus, for any node r, the set 
of points r -I- Ci are also lattice nodes. The divergence at 
a node r is then written as a sum of fluxes deflned on 
the midpoint r + ^c^ of the link connecting the node to 
its neighbour r + . This is schematically represented in 
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Fig. [2] for D2Q9. . Then, Eq. [8] can be discretised as 



E 



(30) 



where Wi are a normahsed set of weights that ensure 
isotropy and and are the deterministic and random 
contributions to the order parameter flux, respectively. 
This ensures the conservation of the order parameter to 
machine accuracy. 

The choice of expressions which relate the fluxes to 
the densities at the nodes must now be dictated by the 
requirement that the FDT holds for all wavevectors. We 
use a symmetric average of node values to compute the 



mid-point fluxes [11 



r = A/i[VAi(r) + V/i(r + c,)] 
-i[(ui^)(r) + (ui^)(r + c,0] 



1 r 
2 



i{r) + i{v + c) 



(31) 
(32) 



Here, c^Sa/s — '^^CiaCij}. To be consistent with this 
choice, the gradient of the chemical potential must be 
computed using 



V/^(r) 



^2 E 



w,Ci^(r + Ci). 



(33) 



parameter. Although this requires, in principle, an al- 
gorithm which updates self-consistently both fields, we 
have to do it sequentially at every time step due to the 
coupling of two different methods, resulting in a hybrid 
scheme for the model H equations. However, we have 
not found any event where the proposed algorithm of al- 
ternate marching in time of FLBE and SMOL leads to 
spurious cross correlations between momentum and order 
parameter fluctuations. 

A number of tests have been carried out to validate the 
algorithm including static and dynamic correlations for 
the order parameter and standard tests for hydrodynam- 
ics. We have always used a D3Q15 model for FLBE with 
lattice units Aa; = At = 1 which leads to a speed of sound 
Cs = To ensure that compressibility is negligible, 

we work in parameter regimes where the Mach number is 
small. Ma ~ u/cg <^ 1. Except when otherwise stated, 
all simulations have been performed on a 32 x 32 x 32 
lattice which is initialized with a uniform random distri- 
bution, and statistics are collected once the system has 
equilibrated. The relaxation parameter, r = 1.1 and 
the temperature kT = 1/3000 are used Q in FLBE un- 
less otherwise specified. Periodic boundary conditions 
are used in all directions in all the simulations. 



A. Order parameter fluctuations 



It is only with the combined choice of the divergence, 
symmetric averaging, and the gradient that the fluctuat- 
ing Cahn-Hilliard equation takes the form 

dMcO + r(q) • (uV')(q) = r(q) . [A/r(q)A*(q)) + |(q)] 

(34) 

where r(q) = WiCi exp(iq-c,;) is the representation of 
the V operator on the lattice. Our choice of discretisa- 
tion ensures that the same operator r(q) appears in both 
the gradient and the divergence in the diffusive term in 
the Cahn-Hilliard equation. As a result, V • V = is 
preserved at all wavevectors, and not only when q — 
as happens with standard discretisations. The resulting 
Laplacian [L(q)]i?v = r(q) • r(q) is less isotropic than 
the Shinozaki-Oono Laplacian as shown in Fig. 15(e) 



in Appendix [Xj Therefore, we use the Shinozaki-Oono 
Laplacian to calculate V^V in the chemical potential. 

Compared to the finite-difference method of the pre- 
vious section, the finite-volume method is not restricted 
to periodic geometries, and thus allows for simulations 
with wall or shear boundary conditions. The compu- 
tational overhead is significantly reduced since the ex- 
pensive Fourier construction of the noise is no longer re- 
quired. 



VI. RESULTS AND VALIDATION 

The order parameter induces a force on the fluid, ac- 
celerating it while the fluid, in turn, advects the order 



We analyze initially a miscible mixture without surface 
tension, characterized by i? = i^T = 0. Since in this case 
the free energy functional, Eq. [2] , is parabolic, the equi- 
librium order parameter distribution follows the Gibbs 
distribution with Gaussian order parameter fluctuations 
of amplitude given in Eq. [181 

Fig. [21 displays the error in the equilibrium ratio (ER) 
between the measured static correlation functions of the 
order parameter and the theoretical prediction, Eq. [TBI 
independent of the wave vector magnitude, as a function 
of the magnitude of the wave vector, for = % = Qz, 
both without and with hydrodynamic coupling. In the 
latter situation we have also compared the performance 
of the finite difference method fscction fY Ap and the finite 
volume method (section [VB[) . In all cases we obtain an 
excellent agreement for the entire wave vector spectrum, 
as opposed to the spurious deviations observed in Fig. 
[H for a standard discretisation of Eq. [51 In Fig. [H the 
velocity-order parameter correlations are plotted using 
both the finite difference and finite volume method to 
show that no spurious scalar-tensor correlations develop 
in the proposed numerical scheme. 

In order to check the homogeneity and isotropy of the 
fiuctuations, polar plots are shown in Fig. [5l In these 
plots, the radius represents the ER as a function of the 
azimuthal angle in a given z plane in lattice space. Dif- 
ferent symbols correspond to three different z planes. ER 
remains essentially unity in all cases, indicating that FDT 
is satisfied in all directions in the lattice. 

The equilibrium structure factor of a miscible binary 
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FIG. 3. (Color Online) Error in the equilibrium ratio as 
a function of wave vector magnitude, q, along the diagonal 
X — y = z considering (1) diffusion alone and (2) coupled 
hydrodynamics with (i) finite difference and (ii) finite volume 
method for the quadratic free energy functional (see Eq. I18|l . 
Simulations have been done on a 32 x 32 x 32 lattice with equi- 
librium initial conditions and parameters used are A = 0.625, 
B = 0, K = 0.0 and M = 0.095. Ensemble averaging is done 
over 10^ time steps and over 25 realizations. 
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FIG. 4. (Color Online) Velocity-order parameter correlation 
for all three components of the velocity in Cartesian coordi- 
nates along the diagonal x = y = z of the domain considering 
coupled hydrodynamics and using the finite difference scheme 
for the quadratic free energy functional (see Eq. I18|l and for 
the same set of parameters as in Fig. [S] Ensemble averaging 
is done over 10* time steps and over 25 realizations. No cross 
correlations are present between fields of different tensorial 
nature. The results obtained using finite volume method are 
shown only for one velocity component for clarity. 



mixture, B — 0, (above the critical temperature) which 
experiences an energy cost to order parameter gradients 
is 



(l^(q)n 



kT 



(35) 



A + Kq^ ' 

On a lattice, the discrete representation of the Laplacian 



must be accounted for, and the static spectrum reads 
accordingly, (|V'(q)P> = kT/{A - KL{q)). 

Since we have used the Shinozaki - Oono form for the 
Laplacian, Eq. IA41 to calculate V'^ip in our simulations, 
— of Eq. [35] is replaced by the Fourier transform of 
appropriate Laplacian I/(q), i.e, Eq. IA8I 

Fig. [6] displays the simulated (|-0(q)p) at equilibrium 
on a wavenumber plane of constant Qz in the absence of 
hydrodynamic coupling while Fig. [7] shows results for 
the full dynamics using the two complementary spatial 
discretisation approaches. The analytical prediction is 
superimposed showing the high degree of accuracy and 
isotropy obtained in all situations. Only at large wave 
vectors the results obtained using the finite difference 
method compare better with theory than those obtained 
from finite volume method. We attribute this accuracy 
loss to the different structure of the lattice Laplacian in 
both approaches, although the errors are consistent with 
the statistical uncertainty associated to the sampling per- 
formed. To show that there is no systematic errors hid- 
den in Fig. [5] and Fig. [3 a one dimensional plot of the 
error in the equilibrium ratio is plotted against q, along 
the diagonal qx = Qy = qz hi the wave vector space in 
Fig. [i 

We have also analyzed the equilibrium dynamic 
structure factor of this miscible mixture, S'(q, t) = 
('(/'(q, t)7/'(q, i + r)), for which we have an analytic ex- 
pression. Taking into account the lattice structure, it 
reads 



kT 



A - KL{q) 



(36) 



Fig. [9] displays In [S'(q, r)/S'(q, 0)] as a function 
of the scaled time T/T(q), where we introduce the 
characteristic decay time for each mode, T(q) ~ 
[M{-L{q)){A - KL{q))]~^. The simulation results re- 
cover the expected slope with a high degree of accuracy 
over all the times covered for each mode for the two 
discretisation schemes of the fiuetuating Cahn-Hilliard 
equation. 



B. Galilean invariance 

The coupling of the order parameter dynamics to the 
fiuid motion must respect Galilean invariance. In order 
to test if the proposed algorithm recovers this basic sym- 
metry, we have imposed a constant velocity along one of 
the system's diagonal, x ~ y. Fig. [TU] displays the or- 
der parameter static structure factor, 5(q) ~ {\tp{q)\'^), 
for a miscible mixture with an energy cost gradient, sub- 
ject to a uniform flow with different magnitudes. Due 
to Galilean invariance, >S'(q) must not be affected by 
the fiuid motion and must coincide with the equilibrium 
curves in Fig. |6] . 

At small flow rates (small Ma), Fig. [TUla, we do not 
see any deviation from the equilibrium predictions, as 
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FIG. 5. (Color Online) Polar plots where the radius indicates the ER as a function of azimuthal angle on lattice points of a 
fixed modulus (16 lattice units from the center), i.e. along a (cos ^, sin 6, z — constant) for (a) diffusion alone, finite difference 
method and (b) hydrodynamics coupled with finite difference method, (c) hydrodynamics coupled with finite volume method. 
Three different symbols o, □, * correspond to 2 = N/8,2N/8 and 3A^/8 planes respectively. Data obtained from the 
same simulations used in Fig. m 




-3-2-10123 




(a) 

FIG. 6. (Color Online) Constant (|V'(g)P) values obtained 
at equilibrium, from a simulation considering diffusion alone 
(without any coupling to hydrodynamics) for the free energy 
functional described by Eq. [2] with B = 0. Simulations ob- 
tained using the finite difference method. Results are shown 
in a wave number plane of (q^ ,qy)- Analytical expression from 
Eq. [33 are superposed onto it using symbols for comparison. 
Simulations are performed on a 32 x 32 x 32 lattice with equi- 
librium initial conditions, A = 0.025, K = O.Of and M = O.f . 
Ensemble averaging is done over f 0^ time steps and over 25 
realizations. 



expected. However, increasing the velocity for Ma > ^, 
Fig. [TUlb shows the development of an anisotropic struc- 
ture factor, which we attribute to the numerical dissi- 
pation associated with advection terms in the order pa- 
rameter conservation equations. Although in principle, 
the proposed LB algorithm does not ensure Galilean in- 
variance at high Ma ( a situation which can be improved 
with complementary LB implementations [39| ) , the main 
source for inaccuracies comes from numerical dissipation 
in the order parameter dynamics. Numerically less dissi- 
pative schemes such as operator splitting may be resorted 
to avoid these limitations (40| . However, in our simu- 
lations we have considered only RK algorithms, which 
recover the correct behavior for small Ma flows. 



C. Fluctuating interfaces 

All the tests described above have used a harmonic 
free energy functional. Below, we present a test of the 
model including the quartic anharmonicity in the free 
energy. At two phase coexistence, with A < and B > 0, 
the order parameter variation across the diffuse interface 
separating the two phases is the well-known hyperbolic 
tangent of Eq. H) In Fig. [TT]we show the order parameter 
profile across the interface, averaged over time and initial 
conditions. We have verified that the mean profile follows 
Eq. 2] with a characteristic width predicted by Eq. \E\ 

Fluctuations about the mean profile are in general 
complicated. However, long-wavelength harmonic fluc- 
tuations are well-described by capillary wave theory 
p3l |4]| . The energy of an interface with instantaneous 
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FIG. 7. (Color Online) Constant values of {\4'{q)\^) from the simulations when the dynamics of the order parameter is coupled 
to the fluid dynamics for the same parameters and lattice size used in Fig. [6] Results for both the finite difference method 
(a) and the finite volume method (b) are shown at a constant plane and expected values from Eq. are superposed as 
symbols. 
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FIG. 8. (Color Online) Error in the equilibrium ratio as 
a function of wave vector magnitude, q, along the diagonal 
qx = qy = qz considering (1) diffusion alone and (2) coupled 
hydrodynamics with (i) finite difference and (ii) finite volume 
method for the free energy functional described by Eq. [2]with 
5 = for the same set of parameters in Fig. |B]and Fig. [T] 



height, h{x,y) is approximated as 




In Fourier space, this is 



from which it follows that 



(IMq)r) = ^. 

7q 



(38) 



(39) 



Since our simulations evolve the entire order parameter 
field, which has both short-wavelength bulk fluctuations 
and long- wavelength capillary fiuetuations, it is necessary 
to tune parameters appropriately to capture the capillary 
fiuetuations. This is ensured when the thermal capillary 

length /* = 



(37) 



4^, the interfacial width and the system 

size A obey P <C Z ^ A. The first inequality ensures 
that the energy scale of the thermal fluctuations excites 
capillary modes and not bulk order parameter modes, 
while the second ensures that the long-wavelength capil- 
lary regime is accessible in the simulation. The capillary 
length condition is equivalently /kT >> 1. 

We have carried out simulations on a system of size 
128 X 1024 where interfaces of linear dimension of 1024 
are symmetrically placed about the center of the domain 
at a gap of 64 lattice units (see Fig. [T2) at two different 
temperatures. Results from these simulations are shown 
in Fig. [T3] using both flnite difference and finite volume 
methods. In diffuse interface models, alternative defini- 



tions of the interface and its location are possible (42| . We 
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-T/x(q) -'C/Tlq) 

(a) (b) 

FIG. 9. (Color Online) Verification of the dynamic correlation function, Eq. 1361 of the order parameter Fourier components. 
Finite difference scheme (a) and finite volume scheme (b) have been used to carry out simulations on a 32 x 32 x 32 lattice 
with A = 0.065, B = 0, K = 0.04 and M = 0.095 with an initial equilibrium distribution. Ensemble averaging is done over 10^ 
time steps and 20 realizations. 




FIG. 10. (Color Online) Galilean invariance of the scheme is tested by applying a uniform velocity field along a diagonal 
direction. Constant values of {|?/'(q)p) from the simulations are plotted along with theoretical predictions as symbols using the 
same parameters as in Fig. [T] (a) At small flow velocities, Ma = 0.08, correct equilibrium is maintained in the simulations, 
(b) However at large flow velocities. Ma = 0.57, an anisotropic distribution of the order parameter fluctuations develops. 



have used a simple linear interpolation to determine the 
location of the interface as the zero of the order param- 
eter. The cross over time for roughening transition and 
the longest relaxation time [i^, Hj] may be estimated as 
~ 10* and ~ 10^ time steps. Therefore, simulation data 
was collected only after 10^ times steps, to ensure sta- 



tionarity of the fluctuations. The logarithmic plot of Fig. 
[T51 shows that the algebraic theoretical prediction can be 
recovered over several orders of magnitude by scaling ap- 
propriately the wave vector and height spectrum magni- 
tudes and changing the system parameters. Exploiting 
the underlying scaling structure of the interface height 
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FIG. 11. (Color Online) Mean equilibrium profile of the or- 
der parameter for two fluid phases coexisting through a planar 
interface. The dashed line is the initial set sharp profile on 
a 64 X 64 lattice of interface in order parameter with left- 
right symmetry (only left half is shown in the plot). Sym- 
bols show the theoretical predictions, (Eq. [4]), the contin- 
uous line is an instantaneous profile from simulations while 
the thick line corresponds to the ensemble averaged profile. 
The continuous line illustrates the magnitude of fluctuations 
around the mean shape. Ensemble averaging is done after 
attaining equilibrium(10^ time steps) over 4 x 10^ time steps 
and 7 realizations. Parameters used in the simulation are 
B ^ 0.025, K = 0.01 and M = 0.1. 
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FIG. 13. (Color Online) Logarithmic plot of the interfacial 
height fluctuation spectrum as a function of the wave vec- 
tor magnitude. Symbols show the simulations results and 
continuous line correspond to the theoretical prediction (Eq. 
I39[) . The wave vector magnitude is scaled with the capil- 
lary length, r , and the magnitude of the height fluctuations 
has also been scaled with to highlight the universal na- 
ture of the capillary spectrum, which is recovered over sev- 
eral orders of magnitude. Four different symbols *, □. x, 
o correspond to simulations with kT = 10"'^ using finite- 
difference method, kT = 10~^ using finite-volume method, 
kT = 1/3000 using finite-difference method and kT = 1/3000 
using finite- volume method respectively, on a 1024 x 128 lat- 
tice (See Fig. I12[l . Free energy and LB simulation parameters 
are -A ^ B ^ 0.05, K = 0.2, M = 0.1 and r = 0.45. 
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FIG. 12. (Color Online) Instantaneous order parameter field 
used for capturing the capillary spectrum. Here two fluid 
phases coexist through two planar fluctuating interfaces. The 
full domain used for simulation is shown on the left side and 
the fluctuating diffused interface is on the right side. The 
continuous line is for -i/i = 0. 



fluctuations, we can combine several numerical simula- 
tions with appropriate fluid parameters to reconstruct 
the whole universal curve, a strategy already exploited in 
the kinetics of phase-separating fluid mixtures |45| . Since 
the quartic anharmonicity is essential in maintaining the 
interface and its fluctuations, this provides a non-linear 



test of the equilibriation in our numerical scheme. 



VII. CONCLUSIONS AND OUTLOOK 

A hybrid method for the numerical solution of the 
model H equations has been developed and validated. A 
fluctuating lattice Boltzmann algorithm is used for hy- 
drodynamics while a stochastic method of lines is pro- 
posed for order parameter conservation equation. Spa- 
tial discretisation in the latter case may be done using 
finite difference or a finite volume schemes both of which 
ensure correct FDT at the lattice level. ELBE takes care 
of fluctuations in momentum at the lattice level. The 
momentum and order parameter equations are coupled 
through stress and advcction terms. The accuracy of 
the algorithm is demonstrated through various hydrody- 
namic and order parameter fluctuation tests. The capil- 
lary spectrum of height fluctuations is reproduced accu- 
rately. 

There are several situations where simulations of fluc- 
tuating hydrodynamics of binary fluid system is neces- 
sary. Eor example, our method can be used to study 
phenomena such as critical fluctuations in symmetric bi- 
nary mixtures and nucleation in asymmetric binary mix- 
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tures. In the light of discussions in section Hill the role of 
noise in the spinodal decomposition of a binary system 
remains unclear [l^. This method may be successfully 
employed in studying the noise driven growth in different 
regimes of the decomposition process. Similarly inter- 
face fluctuations play an important role in several meso 
scale phenomena such as fluctuations driven spreading 



of nano droplets on solid surfaces |12| , d ewetting of thin 
films [31 and break up of nano jets [l3|- Traditionally, 
molecular dynamics simulations have been used to study 
these problems. We expect our mesoscale algorithm to 
be an effective complement to MD simulations which are 
currently limited to short time scales. 
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Appendix A: Comparison of discrete Laplacian 
operators 



In order to ensure the isotropy of the discrete form 
of the Laplacian operator, we have compared four avail- 
able expressions of the operator existing in the literature. 
Controlling the degree of isotropy of the diffusive term of 
the order parameter governing equation (Eq. i.e, V^^ 
is crucial to avoid spurious interface pinning. Note that 
the evaluation of chemical potential (fi) itself contains 
Laplacian of order parameter. We give the details of the 
comparisons here. 

Consider a 3d cubic lattice as shown in Fig. [TH it 
has 6 nearest neighbors, denoted as A^i, 12 next near- 
est neighbors, denoted as N2 and 8 next next nearest 
neighbors, denoted as N3. Correspondingly the set of 
lattice vectors with one, two and three non zero compo- 
nents form the set , c^^ and respectively where 
[cfScf%cf3]gc.. Then 



N3 




FIG. 14. Stencil used for Laplacian calculation for various 
schemes illustrated in this appendix [X] Here TVi is for the 
nearest neighbors, N2 is for next nearest neighbors and N3 is 
for next next nearest neighbors. For clarity only one pair of 
each of them is marked. 



corresponding Fourier transforms are 



12 



30 



30 



6 



12 



22 



i=l 
26 



EV'(r + Ci) - 26?/'(r) 



(Al) 



1 ^ 

-E^(r + c^)-^V'(r) (A3) 



(A4) 



where V^r) = i^{x,y,z). The sufSxes CD, PK, SO 
and LB stand for central difference, Patra-Kartunnen, 
Shinozaki-Oono and lattice Boltzmann, respectively. Eq. 
IA2I is the standard central finite difference expression. 
Eq. IA3I has been systematically derived by imposing 
conditions of rotational invariance and isotropy of the 
operator [4^ . Eq. IA4I is popular in the cell-dynamics 
and phase separation studies |3a |. Eq. IA4l is a simple ex- 
pression used in lattice Boltzmann simulations (47[ . The 



c.]-3} 



(A5) 



1 



[i(q)]pK = ^{28 [c. 



+12 [c^cy + c^c, + Cyc] + 8 [c:,cyc^] " 128}(A6) 



[^(q)]so = ^{12 [c, + c,y + c,] 



+12 [cxCy + Cj,Cz + CyCz] + 8 [cxCyC^] ~ 80} (A7) 
1 
9 

+4[c xCy ~\~ CxCz H~ CyCzl H~ 8 \CxCyCz ] - 26}(A8) 



[H<i)]LB = ^ {2 [c^ + + C;, 



where = cos 



cos 



Clearly all the Laplacian operators are negative defi- 
nite (except at g = 0). Fig. [12] shows the magnitude 
of the different expressions of the Laplacian operator in 
wave vector planes with constant Qz- These plots clearly 
display the four fold symmetry of the lattice. Nonethe- 
less, the effect is less pronounced for th e expression sug- 
gested by Shinozaki and Oono (15(c)) and so we have 
used Eq. IA4I for the calculations in section |V] - 1 VII 

In the finite volume approach to solve the order param- 
eter evolution, fluxes are calculated (Eq. on the links 
connecting the lattice nodes. Ensuring FDT leads to an 
equivalent Laplacian operator whose Fourier transform 
reads 



1 r 2 
[i(q)]Fy = -g |[2s:i. + SxCyC^] +[2 

+ [25^ + SzCa:Cyf^ 



Sy -T SyCxCz] 



(A9) 



where Sx = sin qx,Sy = sin qy,Sz = sing^. This is also 
plotted in Fig. [12] for comparison purpose. 
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Appendix B: Fluctuating Navier-Stokes solver 

In this appendix we describe FLBE method we use 
to solve the fluctuating Navicr-Stokcs equations. A self- 
contained derivation is given below to include noise 
and external force densities [l^l to the standard LB 
model. 

In a standard DdQn LEE model where the velocity is 
discretized into n components in a d dimensional space, 
the discrete form of the fluctuating Boltzmann equation 
is given in Eq. 1191 where the moments of the single par- 
ticle distribution function fi, are expressed in Eq. 1201 
A multi-scale expansion, or a moment closure method, 
shows that the above equation has Eq. [9] as its hydro- 
dynamic limit |3l| . Since FLBE is a hyperbolic equa- 



tion with local non-linearities, it is considerably easier to 
solve than Eq. [HI which has a parabolic-hyperbolic char- 
acter with advective non-linearities. The methodology of 
the FLBE has been explained in detail in Q, while the 
method by which force densities are added is given in de- 
tail in [s^l • Here we outline the integration scheme we use 
when force densities and fluctuating forces are combined 
in the FLBE. 

We can rearrange equation, Eq. [TO] to obtain 

at/^ + c,.V/, = i?,(x,t) (Bl) 

where Ri{x,t) = -J^j^ijifj ~ fj) + represents 
the effects of collision, forcing and thermal fluctuations. 



<I>i = Q — F ■ V cf accounts for the fluctuating and ex- 
ternal forces acting on the distribution function. Using 
the method of characteristics, this set of first order hy- 
perbolic equations can be integrated over a time interval 
A< to get 

/i(x + CiAt,i + At)-/i(x,i) = / dsi?i(x + CiS,f + s) 

Jo 

(B2) 

The integral above may be approximated to second order 
accuracy using the trapezium rule and the resulting terms 
transposed to give a set of implicit equations for the fi : 

f,{yi + CiAt,t + At) - ^i?,(x + c,Ai,i + At) = 

/,(x,t)-^i?,(x,t)+Ati?,(x,t). (B3) 

In terms of the auxiliary distribution function, Eq. I2H 
the evolution equation reduces to 

/,(x-fc,Ai,t + At) ==/,(x,t)+i?,(x,t)At. (B4) 

indicating that we can understand LBE evolution 
through a simple relaxational step in which the dis- 
tributions fi are relaxed to their postcoUisional values 

/»(x,r), 



/,(x,r) = /,(x,t)+i?,(x,t)At, 



(B5) 
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followed by a propagation step in which the postcolli- 
sional distributions are propagated along a Lagrangian 
trajectory without further change, 



/,(x + c,At,t + At) ^ /»(x,r). 



(B6) 



Thus the computational part of the method is most nat- 
urally framed in terms of the auxiliary distributions fi 
instead of the physical distribution functions fi them- 
selves. To obtain the postcollisional without having to 
refer to the ft, the latter must be eliminated from Eq. 
IB4I Inverting the equations defining the fi in Eq. [21] we 
obtain 



= ( 1 + [-LMfk - fk) + '^A^,t)]- (B7) 



Combining this with Eq. IB4I we obtain a numerical 
scheme for the discrete Boltzmann equation with a gen- 
eral collision operator in terms of the fi'. 



/,(x + c, At, t + At) = /,(x, t) + 



At 



1 + [-L,fc(A.-/l?) + <i>,(x,t)]At. (B8) 



For a single time relaxation operator, where Lij ~ Sij /t, 



this takes on a particularly simple form, 
/,(x + c,At,t + At) = /,(x,t) 



At 



r + At/2 



h(/,-/°) + r$.(x,t)], 



(B9) 



For a nondiagonal collision operator, the collision term is 
best evaluated in the moment basis. For example, using a 
collision operator in which the ghost modes are projected 
out [sij and the stress modes relax at a rate of t~^, the 
post coUisional fi are given by 



/i(x,t*) ^ Wt[ p + 



2c4 



(BIO) 



where the normalized weights Wi ensure the isotropy, and 
Aa , the momentum component of the postcollisional aux- 
iliary distributions, is 



i=0 



pFaAt 



(Bll) 



while Bai3, the stress component, reads 

At 



Bal3 — ^ fiQ 



ial3 



^ ^ fiQiaP 

4=0 



T + At/2 

+pVaVl3 + riVaFp + FaVp) + T ^ CiQiafl -(612) 



i=0 



The mass and momentum densities are obtained as p = 
SLo/i pva = J2i=ofiCia + pFa^, respectively. 
The equilibria can be reconstructed from p and pv. 



